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Abstract. Given a "black box" function to evaluate an unknown ra- 
tional polynomial f e <Q[x] at points modulo a prime p, we exhibit algo- 
rithms to compute the representation of the polynomial in the sparsest 
shifted power basis. That is, we determine the sparsity t G Z>0) the 
shift CK G Q, the exponents < ei < 62 < ■ ■ ■ < ej, and the coefHcients 
ci , . . . , ct G Q \ {0} such that 

f{x) = ci{x - +C2{x - -\ \- ct{x - ay*. 

The computed sparsity t is absolutely minimal over any shifted power 
basis. The novelty of our algorithm is that the complexity is polyno- 
mial in the (sparse) representation size, which may be logarithmic in 
the degree of /. Our method combines previous celebrated results on 
sparse interpolation and computing sparsest shifts, and provides a way 
to handle polynomials with extremely high degree which are, in some 
sense, sparse in information. 
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1. Introduction 

Interpolating an unknown polynomial from a set of evaluations is a problem 
which has interested mathematicians for hundreds of years, and which is now 
implemented as a standard function in most computer algebra systems. To 
illustrate some different kinds of interpolation problems, consider the following 
three representations for a polynomial / of degree n: 

(1.1) f{x) = ao + aix + a2x'^ -\ h a„x", 

(1.2) f(x) = bo + hx"^' + b2X'^' + --- + bsX'^% 

(1.3) f{x) = Co + ci{x - ay + C2{x - ay^ ^ hct{x-ay*. 
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In (1.1) we see the dense representation of the polynomial, where all coefficients 
(even zeroes) are represented. Newton and Waring discovered methods to in- 
terpolate / in time proportional to the size of this representation in the 18*^^ 
century. The sparse or lacunary representation is shown in (1.2), wherein only 
the terms with non-zero coefficients are written (with the possible exception 
of the constant coefficient bo). Here we say that / is s- sparse because it has 
exactly s non-zero and non-constant terms; the constant coefficient requires 
special treatment in our algorithms regardless of whether or not it is zero, and 
so we do not count it towards the total number of terms. Ben-Or & Tiwari 
(1988) discovered a method to interpolate in time polynomial in the size of this 
representation. Kaltofen & Lee (2003) present and analyze very efficient algo- 
rithms for this problem, in theory and practice. The Ben-Or & Tiwari method 
has also been examined in the context of approximate (floating point) polyno- 
mials by Giesbrecht et al. (2006), where the similarity to the 1795 method of de 
Prony is also pointed out. Blaser et al. (2009) consider the more basic problem 
of identity testing of sparse polynomials over Q, and present a deterministic 
polynomial-time algorithm. 

In (1.3), / is written in the shifted power basis 1, (x — a), (x — a)^, . . ., and 
we say that a is a t-sparse shift of / because the representation has exactly t 
non-zero and non-constant terms in this basis. When a is chosen so that t is 
absolutely minimal in (1.3), we call this the sparsest shift of /. We present new 
algorithms to interpolate / G Q[x], given a black box for evaluation, in time 
proportional to the size of the shifted-lacunary representation corresponding 
to (1.3). It is easy to see that t could be exponentially smaller than both n 
and s, for example when f = [x + 1)", demonstrating that our algorithms are 
providing a significant improvement in complexity over those previously known, 
whose running times are polynomial in n and s. 

The main applications of all these methods for polynomial interpolation 
are signal processing and reducing intermediate expression swell. Dense and 
sparse interpolation have been applied successfully to both these ends, and 
our new algorithms effectively extend the class of polynomials for which such 
applications can be made. 

The most significant challenge here is computing the sparsest shift a G Q. 
Computing this value from a set of evaluation points was stated as an open 
problem by Borodin & Tiwari (1991). An algorithm for a generalization of 
our problem in the dense representation was given by Grigoriev & Karpinski 
(1993), though its cost is exponential in the size of the output; they admit that 
the dependency on the degree of the polynomial is probably not optimal. Our 
algorithm achieves deterministic polynomial-time complexity for polynomials 
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over the rational numbers. We are always careful to count the bit complexity 
— the number of fixed-precision machine operations — and hence account for 
any coefficient growth in the solution or intermediate expressions. 

The black box model we use is slightly modified from the traditional one: 

pen,e eZp — — ^ f{e) mod p 

fix) G Q[x] 

Given a prime p and an element 6 in Zp, the black box computes the value of 
the unknown polynomial evaluated at 6 over the field Zp. (An error is produced 
exactly in those unfortunate circumstances that p divides the denominator of 
f{9).) We generally refer to this as a modular black box. To account for the 
reasonable possibility that the cost of black box calls depends on the size of 
p, we define kj to be an upper bound on the number of field operations in Zp 
used in black box evaluation, for a given polynomial / G Q[a;]. 

Some kind of extension to the standard black box, such as the modular 
black box proposed here, is in fact necessary, since the value of a polynomial 
of degree n at any point other than 0, ±1 will typically have n bits or more. 
Thus, any algorithm whose complexity is proportional to log n cannot perform 
such an evaluation over Q or Z. Other possibilities might include allowing 
for evaluations on the unit circle in some representation of a subfield of C, or 
returning only a limited number of bits of precision for an evaluation. 

To be precise about our notion of size, first define size(g) for g G Q to be 
the number of bits needed to represent q. So if we write q = f with a G Z, 
6 G N, and gcd(a, b) = 1, then size(g) = [log2(|a| + 1)] + [log2(6 + 1)] + 1. For 
a rational polynomial / as in (1.3), define: 

t t 

(1.4) size(/) = size(Q;) + ^^size(ci) + ^^size(ej). 

i=0 i=l 

We will often employ the following upper bound for simplicity: 

(1.5) size(/) < size(a) + t {H{f) + log, n) , 

where H{f) is defined as maxo<i<t size(ci). 

Our algorithms will have polynomial complexity in the smallest possible 
size(/). For the complexity analysis, we use the normal notion of a "mul- 
tiplication time" function M{n), which is the number of field operations re- 
quired to compute the product of polynomials with degrees less than n, or 
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integers with sizes at most n. We always assume that M(n) G ^7(n) and 
M(n) G O(ra^). Using the results from Cantor & Kaltofen (1991), we can 
write M(n) G 0(nlognloglog?T,). 

The remainder of the paper is structured as follows. In Section 2 we show 
how to find the sparsest shift from evaluation points in Zp, where p is a prime 
with some special properties provided by some "oracle" . In Section 3 we show 
how to perform sparse interpolation given a modular black box for a poly- 
nomial. In Section 4 we show how to generate primes such that a sufficient 
number satisfy the conditions of our oracle. Section 5 provides the complexity 
analysis of our algorithms. We conclude in Section 6, and introduce some open 
questions. 

2. Computing the Sparsest Shift 

For a polynomial / G Q[x], we first focus on computing the sparsest shift a G Q 
so that f{x + a) has a minimal number of non-zero and non-constant terms. 
This information will later be used to recover a representation of the unknown 
polynomial. 

2.1. The polynomial f^'^\ Here, and for the remainder of this paper, for 
a prime p and / G Q[a;], define /(p) G Zp[s] to be the unique polynomial 
with degree less than p which is equivalent to / modulo — x and with all 
coefficients reduced modulo p. From Fermat's Little Theorem, we then see 
immediately that f^P\a) = /(a)modp for all a G Zp. Hence f^P^ can be 
found by evaluating / at each point 0, 1, ... ,p — 1 modulo p and using dense 
interpolation over Zp[x]. 

Notice that, over Zp[x], {x — aY = x — a mod x^ — x, and therefore {x — 
aY' = [x — aY for any k ^ such that Ci = k mod {p — 1). The smallest such 
k is in the range {1,2, . . . we now define this with some more notation. 
For a G Z and positive integer m, define a remi m to be the unique integer in 
the range {1, 2, . . . , m} which is congruent to a modulo m. As usual, aremm 
denotes the unique congruent integer in the range {0, 1, . . . , m — 1}. 

If / is as in (1.3), then by reducing term- by-term we can write 

t 

(2.1) fP\x) = (coremp) + ^(ciremp)(x - apY'''''^'^P-^\ 

1=1 

where Up is defined as a rem p. Hence, for some k < t, ap is a /c-sparse shift 
for f^P\ That is, the polynomial f^\x + ap) over Zp[x] has at most t non-zero 
and non-constant terms. 
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Computing f^^^ from a modular black box for / is straightforward. First, 
use p black-box calls to determine f{i) remp for i = 0,1, . . . ,p — 1. Recalling 
that is the number of field operations in Zp for each black-box call, the cost 
of this step is 0(pK/M(logp)) bit operations. Second, we use the well-known 
divide-and-conquer method to interpolate f^^^ into the dense representation 
(see, e.g., Borodin & Mimro (1975, Section 4.5)). Since degf^^^ < p, this step 
has bit complexity 0(M(p)M(logp) logp). 

Furthermore, for any a G Zp, the dense representation of f^^\x + a) can be 
computed in exactly the same way as the second step above, simply by shifting 
the indices of the already-evaluated points by a. This immediately gives a 
naive algorithm for computing the sparsest shift of f^^^: compute f^^\x + 'y) for 
7 = 0, 1, . . . ,p—l, and return the 7 that minimizes the number of non-zero, non- 
constant terms. The bit complexity of this approach is 0{p\ogp M(p)M(logp)), 
which for our applications will often be less costly than the more sophisticated 
approaches of, e.g., Lakshman & Saunders (1996) or Giesbrecht et al. (2003), 
precisely because p will not be very much larger than deg/^^\ 

2.2. Overview of Approach. We will make repeated use of the following 
fundamental theorem from Lakshman & Saunders (1996): 

Fact 2.2. Let F be an arbitrary field and f G F[x], and suppose a E F is such 
that f{x + a) has t non-zero and non-constant terms. If deg / > 2t + 1 then a 
is the unique sparsest shift of f . 

From this we can see that, if a is the unique sparsest shift of /, then 
ap = aremp is the unique sparsest shift of f^^^ provided that deg/^^^ > 2t-|- 1. 
This observation provides the basis for our algorithm. 

The input to the algorithms will be a modular black box for evaluating a 
rational polynomial, as described above, and bounds on the maximal size of 
the unknown polynomial. Note that such bounds are a necessity in any type 
of black-box interpolation algorithm, since otherwise we could never be sure 
that the computed polynomial is really equal to the black-box function at every 
point. Specifically, we require Ba, Bt, Bh, Bn G N such that 

size(Q;) < Ba, 

t<BT, 

size(cj) < Bh, for < i < t, 
\0g2n < Bn- 
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By considering the following polynomial: 

c{x - a)" + {x- + ■ ■ ■ + (x - 

we see that these bounds are independent — that is, none is polynomially- 
bounded by the others — and therefore are all necessary. 

We are now ready to present the algorithm for computing the sparsest shift 
a almost in its entirety. The only part of the algorithm left unspecified is 
an oracle which, based on the values of the bounds, produces primes to use. 
We want primes p such that deg/*^^^ > 2t + 1, which allows us to recover one 
modular image of the sparsest shift a. But since we do not know the exact value 
of t or the degree n of f over Q[x], we define some prime p to be a good prime 
for sparsest shift computation if and only if deg /^^^ > min{2i?r + l, n}. For the 
remainder of this section, "good prime" means "good prime for sparsest shift 
computation." Our oracle indicates when enough primes have been produced 
so that at least one of them is guaranteed to have been a good prime, which is 
necessary for the procedure to terminate. The details of how to construct such 
an oracle will be considered in Section 4. 

Algorithm 2.3. Computing the sparsest shift. 

Input: o A modular black box for an unknown polynomial / G Q[a;] 

o Bounds Ba, Bt, Bh, -Bat G N as described above 

o An oracle which produces primes and indicates when at least one 
good prime must have been produced 

Output: A sparsest shift a of /. 

1. 1, ^ ^0 

2. While log2 P < 2Ba + 1 do 3-14 

3. p ^ new prime from the oracle 

4. Evaluate f{i) remp for i = 0,1, . . . ,p — 1 

5. Use dense interpolation to compute f^^^ 

6. If deg /(P) > 2Bt + 1 then 

7. Use dense interpolation to compute f^^\x + 'j) for 7 = 1, 2, . . . ,p— 1 

8. ttp ^ the unique sparsest shift of f^^^ 

9. p^p-p, g<-g[j{p} 

10. Else if P = 1 and oracle indicates > 1 good prime has been produced 
then 

11. q ^ least prime such that log2 q > 2BtBa+Bh (computed directly) 

12. Evaluate f{i) rem g for i = 0, 1, . . . , 2Bt 
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13. Compute f E Q[x] with deg / < 2Bt by dense interpolation in Zq[x] 
followed by rational reconstruction on the coefficients 

14. Return A sparsest shift a computed by a univariate algorithm from 
Gicsbrccht et al. (2003) on input / 

15. Return The unique a = a/b e Q such that \a\,b < 2^^ and a = bap mod p 
for each p E Q, using Chinese remaindering and rational reconstruction 

Theorem 2.4. With inputs as specified, Algorithm 2.3 correctly returns a 
sparsest shift a of f . 

Proof. Let /, Ba, Bt, Bh, Bm be the inputs to the algorithm, and suppose 
t,a are as specified in (1.3). 

First, consider the degenerate case where n < 25^, i.e., the bound on the 
sparsity of the sparsest shift is at least half the actual degree of /. Then, 
smce each /(P) can have degree at most n (regardless of the choice of p), the 
condition of Step 6 will never be true. Hence Steps 10-14 will eventually be 
executed. The size of coefficients over the standard power basis is bounded by 
2BtBa + Bh since deg / < 2Bt, and therefore / will be correctly computed 
on Step 5. In this case. Fact 2.2 may not apply, i.e. the sparsest shift may not 
be unique, but the algorithms from Giesbrecht et al. (2003) will still produce 
a sparsest shift of /. 

Now suppose instead that n > 2Bt + 1. The oracle eventually produces a 
good prime p, so that deg f^^^ > 2Bt + 1. Since t < Bt and /(^^ has at most t 
non-zero and non-constant terms in the (aremp)-shifted power basis, the value 
computed as ctp on Step 8 is exactly aremp, by Fact 2.2. The value of P will 
also be set to p > 1 here, and can only increase. So the condition of Step 10 
is never true. Since the numerator and denominator of a are both bounded 
above by 2^^., we can use rational reconstruction to compute a once we have 
the image modulo P for P > 2^-^^+^. Therefore, when we reach Step 15, we 
have enough images ap to recover and return the correct value of a. □ 

We still need to specify which algorithm to use to compute the sparsest 
shift of a densely-represented / e Q[x] on Step 14. To make Algorithm 2.3 
completely deterministic, we should use the univariate symbolic algorithm from 
Giesbrecht et al. (2003, Section 3.1), although this will have very high complex- 
ity. Using a probabilistic algorithm instead gives the following, which follows 
directly from the referenced work. 
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Theorem 2.5. If the "two projections" algorithm of Giesbrecht et al. (2003, 
Section 3.3) is used on Step 14, then Steps 10-14 of Algorithm 2.3 can be per- 
formed with 0{B^M{B^Ba + B^Bh)) bit operations, plus 0{KfBTM{BTBA + 
Bh)) bit operations for the black-box evaluations. 

The precise complexity analysis proving that the entire Algorithm 2.3 has 
bit complexity polynomial in the bounds given depends heavily on the size and 
number of primes p that are used, and so must be postponed until Section 5.1, 
after our discussion on choosing primes. 

Example 2.6. Suppose we are given a modular black box for the following 
unknown polynomial: 

f(x) = x^^ - 45x^^ + 945x^^ - 12285x^^ + 110565a;^^ - 729729x^° 

+ 3648645a;^ - 14073345a;^ + 42220035a;^ - 98513415x^+ 
177324145x^ - 241805625x^ + 241805475x^ - 167403375x2 
+ 71743725a; - 14348421, 

along with the bounds Ba = 4, Brp = 2, B^ = 4, and B^ = 4. One may easily 
confirm that f{x) = (x — 3)^^ — 2 (x — 3)^, and hence these bounds are actually 
tight. 

Now suppose the oracle produces p = 7 in Step 3. We use the black box to 
find /(O), /(I), . . . , /(6) in Z7, and dense interpolation to compute 

= 5x^ + 2x^ + 3x^ + 6^^ + x + 4. 

Since deg/^''^ = 5 > 2Bt + 1, we move on to Step 8 and compute each 
f^'^\x + 7) with 7 = 1, 2, . . . , 6. Examining these, we see that /'•^^(x + 3) = 
5x^ + x^ has the fewest non-zero and non-constant terms, and so set ay to 3 
on Step 8. This means the sparsest shift must be congruent to 3 modulo 7. 
This provides a single modular image for use in Chinese remaindering and ra- 
tional reconstruction on Step 15, after enough successful iterations for different 
primes p. O 

2.3. Conditions for Success. We have seen that, provided deg/ > 2Bt, a 
good prime p is one such that deg /^^^ > 2Bt. The following theorem provides 
(quite loose) sufficient conditions on p to satisfy this requirement. 

Theorem 2.7. Let / G Q[x] as in (1.3) and Bt eN such that t < Bt. Then, 
for some prime p, the degree of f^^^ is greater than 2Bt whenever the following 
hold: 
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o Q ^ mod p; 

o Vi e {1, 2, . . . , t - 1}, et ^ Cj mod {p - 1); 
o Vi G {!,..., 25t}, et ^ i mod (p- 1). 

Proof. The first condition guarantees that the last term of f^^\x) as in 
(2.1) does not vanish. We also know that there is no other term with the same 
degree from the second condition. Finally, the third condition tells us that the 
degree of the last term will be greater than 2Bt- Hence the degree of /^^^ is 
greater than 2Bt- □ 

For purposes of computation it will be convenient to simplify the above 
conditions to two non-divisibility requirements, on p and p — 1 respectively: 

Corollary 2.8. Let /, Bt, Bh, Bn he as in the input to Algorithm 2.3 with 
deg/ > 2Bt. Then there exist Ci, C2 G N with log2 Ci < 2Bh and logs C2 < 
Bn{Wt - 1) such that deg/(P) > 2Bt whenever p \ Ci and (p - 1) f C2. 

Proof. Write / as in (1.3). We will use the sufficient conditions given in 
Theorem 2.7. Write \ct\ = a/b for a, 6 G N relatively prime. In order for q remp 
to be well-defined and not zero, neither a nor b can vanish modulo p. This is 
true whenever p f ab. Set Ci = ab. Since a,b < 2^" , logs C*! = log2(a6) < 2Bh. 
Now write 

t-l 2Bt 

C2 = l[{et-ei)-l[{et-i). 

1=1 1=1 

We can see that the second and third conditions of Theorem 2.7 are satisfied 
whenever (p — 1) f C2. Now, since each integer Cj is distinct and positive, and Ct 
is the greatest of these, each (e^ — Cj) is a positive integer less than Cf. Similarly, 
since Ct = deg / > 2i?y, each (ct — i) in the second product is also a positive 
integer less than Therefore, using the fact that t < Bt, we see C2 < e\^'^~^ . 
Furthermore, < 2^^, so we know that logjCs < Bn{3Bt — 1). □ 

A similar criteria for success is required in Blaser et al. (2009), and they 
employ Linnik's theorem to obtain a polynomial-time algorithm for polynomial 
identity testing. Linnik's theorem was also employed in Gicsbrecht & Roche 
(2007) to yield a much more expensive deterministic polynomial-time algorithm 
for finding sparse shifts than the one presented here. 
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3. Interpolation 

Once we know the value of the sparsest shift a of /, we can trivially construct 
a modular black box for the t-sparse polynomial f{x + a) using the modular 
black box for /. Therefore, for the purposes of interpolation, we can assume 
a = 0, and focus only on interpolating a t-sparse polynomial / e Q[x] given a 
modular black box for its evaluation. The basic techniques of this section are, 
for the most part, known in the literature. However, a unified presentation in 
terms of bit complexity for our model of modular black boxes will be helpful. 
For convenience, we restate the notation for / and f^P\ given a prime p: 

(3.1) / = Co + cix^' + C2X^' + ■■■ + CtX'\ 

(3.2) /(P) = (coremp) + {ciTemp)x^'''''^'^P-^^ + ■■■ + {ct rem p)x^''''''''^f-^\ 

Again, we assume that we are given bounds Bh, Bt, and i?7v on maxj size(cj), 
t, and log2deg/, respectively. We also introduce the notation r(/), which is 
defined to be the number of distinct non-zero, non-constant terms in the uni- 
variate polynomial /. 

This algorithm will again use the polynomials /*^^^ for primes p, but now 
rather than a degree condition, we need f^^^ to have the maximal number of 
non-constant terms. So we define a prime p to be a good prime for interpolation 
if and only if r^f^^^) = t. Again, the term "good prime" refers to this kind of 
prime for the remainder of this section. 

Now suppose we have used modular evaluation and dense interpolation (as 
in Algorithm 2.3) to recover the polynomials Z^^-* for k distinct good primes 
Pi,...,Pk. We therefore have k images of each exponent Cj modulo {pi — 
1) , . . . , {pk — 1) . Write each of these polynomials as: 

(3.3) /(^•)=c« + cSV%--- + clV''. 

Note that it is not generally the case that e^*^ = ejTemi{pi — 1). Because we 
don't know how to associate the exponents in each polynomial Z^^'-* with their 
pre-image in Z, a simple Chinese remaindering on the exponents will not work. 
Possible approaches are provided by Kaltofcn (1988), Kaltofen et al. (1990) or 
Avendaiio et al. (2006). However, the most suitable approach for our purposes 
is the clever technique of Garg & Schost (2009), based on ideas of Grigoriev & 
Karpinski (1987). We interpolate the polynomial 

(3.4) g{z) = {z- ei){z - 62) ■ ■ ■ (z - Ct), 

whose coefficients are symmetric functions in the Cj's. Given f^^'\ we have all 
the values of e^*'* remi (pj — 1) for j = 1, . . . ,t; we just don't know the order. 
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But since g is not dependent on the order, we can compute g mod (pj — 1) for 
i = 1, . . . ,k, and then find the roots of g E to determine the exponents 
ei, . . . ,et. Once we know the exponents, we recover the coefficients from their 
images modulo each prime. The correct coefficient in each f^^^ can be identified 
because the residues of the exponents modulo p — 1 are unique, for each chosen 
prime p. This approach is made explicit in the following algorithm. 

Algorithm 3.5. Sparse Polynomial Interpolation over Q[x]. 
Input: o A modular black box for unknown / G Q[a;] 
o Bounds Bh and Bj^i as described above 

o An oracle which produces primes and indicates when at least one 
good prime must have been returned 

Output: / G Q[x] as in (3.1) 

1. g ^ 1, P ^ 1, k^l, t^O 

2. While log2 P < 2Bh + 1 or loga Q < Bn 

or the oracle does not guarantee a good prime has been produced 
do 3-8 

3. Pk ^ new prime from the oracle 

4. Compute /(p*:) by black box calls and dense interpolation 

5. If r(/(P^)) > t then 

6. Q^pk-l,P^ Pk, t i- r(/(f'=)), p, ^ pk, /(P^) <- k ^ 2 

7. Else if T{fP^^) = t then 

8. Q ^ lcm((5,pfc - 1), P^P-pk, k^k + 1 

9. For i G {1,...,A;- 1} do 

10. gip.) ^ Yl^^.^^(z - ef) mod Pi -1 

11. Construct g = aQ+aiZ+a2Z^+- ■ ■+atZ^ G Z[x] such that g = f?'-^'-' mod Pi — 1 
for 1 < i < A;, by Chinese remaindering 

12. Factor g as {z — ei){z — 62) ■ ■ ■ {z — to determine ei, . . . , Ct G Z 

13. For 1 < i < t do 

14. For 1 < j < A; do 

15. Find the exponent ef^ of /^^^^ such that e'^^ = Cj mod pj — 1 

16. Reconstruct Cj G Q by Chinese remaindering from residues , . . . , c^'^J 

17. Reconstruct cq G Q by Chinese remaindering from residues c^q \ . . . , c^^ 

The following theorem follows from the above discussion. 
Theorem 3.6. Algorithm 3.5 works correctly as stated. 
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Again, this algorithm runs in polynomial time in the bounds given, but we 
postpone the detailed complexity analysis until Section 5.2, after we discuss 
how to choose primes from the "oracle". Some small practical improvements 
may be gained if we use Algorithm 3.5 to interpolate f{x + a) after running 
Algorithm 2.3 to determine the sparsest shift a, since in this case we will have 
a few previously-computed polynomials f^P\ However, we do not explicitly 
consider this savings in our analysis, as there is not necessarily any asymptotic 
gain. 

Now we just need to analyze the conditions for primes p to be good. This is 
quite similar to the analysis of the sparsest shift algorithm above, so we omit 
many of the details here. 

Theorem 3.7. Let f,BT,BH,BN be as above. There exist Ci,C2 G N with 
logaCi < 2BhBt and log2C2 < \BnBt{Bt - 1) such that T{f p'^) is maximal 
whenever p\ Ci and {p — 1) \ €2- 

Proof. Let / be as in (3.1), write \ci\ = ai/bi in lowest terms for i = 1, . . . ,t, 
and define 

t t t 

Ci = l[aA, C2 = l[l[ (e,-e,). 

i=l 1=1 j=i+l 

Now suppose p is a prime such that p \ Ci and (p — 1) f C2. From the 
first condition, we see that each q mod p is well-defined and nonzero, and so 
none of the terms of /'•^•' vanish. Furthermore, from the second condition, 
Ci ^ Cfc mod p — 1 for all i 7^ j, so that none of the terms of collide. 
Therefore /'•^•' contains exactly t non-constant terms. The bounds on Ci and 
C2 follow from the facts that each a^, 6j < 2^" and each difference of exponents 
is at most 2^^. □ 

4. Generating primes 

We now turn our attention to the problem of generating primes for the sparsest 
shift and interpolation algorithms. In previous sections we assumed we had an 
"oracle" for this, but now we present an explicit and analyzed algorithm. 

The definition of a "good prime" is not the same for the algorithms in 
Section 2 and Section 3. However, Corollary 2.8 and Theorem 3.7 provide a 
unified presentation of sufficient conditions for primes being "good". Here we 
call a prime which satisfies those sufficient conditions a useful prime. So every 
useful prime is good (with the bounds appropriately specified for the relevant 
algorithm), but some good primes might not be useful. 
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We first describe a set V of primes sucli tliat tlie number and density of 
useful primes within tlie set is sufficiently higli. We will assume tliat tliere exist 
numbers Ci, C2, and useful primes p are tfiose sucli tliat p\ Ci and (p — 1) f C2. 
The numbers Ci and C2 will be unknown, but we will assume we are given 
bounds f3i, (32 such that logg Ci < /3i and log2 C2 < (32- Suppose we want to 
find i useful primes. We construct V explicitly, of a size guaranteed to contain 
enough useful primes, then enumerate it. 

The following fact is immediate from Mikawa (2001), though it has been 
somewhat simplified here, and the use of (unknown) constants is made more 
explicit. This will be important in our computational methods. 

For g G Z, let S{q) be the smallest prime p such that q\{p — 1). 

Fact 4.1 (Mikawa 2001). There exists a constant n > 0, such that for all n > 
H, and for all integers q G {n, . . . , 2n} with fewer than /in/ log^ n exceptions, 
we have S{q) < q^-^^. 

Our algorithms for generating useful primes require explicit knowledge of 
the value of the constant /i in order to run correctly. So we will assume that 
we know /i in what follows. To get around the fact that we do not, we simply 
start by assuming that /i = 1, and run any algorithm depending upon it. If 
the algorithm fails we simply double our estimate for /i and repeat. At most a 
constant number of doublings is required. We make no claim this is particularly 
practical. 

For convenience we define 

T(x) = - 

5 log X log^ X 

Theorem 4.2. Let logaCi < log2C2 < /32 and £ be as above. Let n be 
the smallest integer such that n > 21, n > fi and T(n) > f3i + (32 + i. Define 

Q = {q prime : n < q < 2n and S{q) < q^'^^}, V = {S{q) : g G Q}. 

Then the number of primes in V is at least (3i + (32 + I, and the number of 
useful primes in V , such that p\C\ and (p — 1) f C2, is at least i. For all p eV 
we have p G 0((/3i + (32 + £)^-«^ ■ log^-'9(/3i + (32 + £)). 

Proof. By Rosser & Schoenfeld (1962), the number of primes between n 
and 2n is at least 3?7,/(51ogn) for n > 21. Applying Fact 4.1, we see #Q > 
3n/(51ogn) — fin/log^n when n > max{yU, 21}. Now suppose ^(gi) = S{q2) 
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for qi,q2 G Q. If qi < q2, then S{qi) > qf, a contradiction with the definition 
of Q. So we must have qi = q2, and hence 

#P = #g>T(n)>/3i + /32 + £. 

We know that there are at most logg Ci < P2 primes p & V such that p\Ci. We 
also know that there are at most log2 C2 < ^2 primes g G Q such that g | C2, 
and hence at most logg C2 primes p eV such that p = S{q) and g | (p — 1) | Ci. 
Thus, by construction V contains at most Pi + (32 primes that are not useful 
out of (3i + (32 + ^ total primes. 

To analyze the size of the primes in V, we note that to make T(ra) > 
Pi + (32 + ^, we have n G &{{(3i + (32 + i) ■ log(/3i + (32 + i)) and each q e Q 
satisfies q G 0{n). Elements of V will be of magnitude at most (2?7,)^ '^^ and 
hence p G 0{{(3i + (32 + i)'-^^ log'-''(/3i + (32 + e)). □ 

Given (3i, (32 and i as above (where logj Ci < (3i and log2 C2 < (32 for 
unknown Ci and C2), we generate the primes in V as follows. 

Start by assuming that /i = 1, and compute n as the smallest integer such 
that T(n) >Pi + (32+i,n>^ and n > 21. List all primes between n 
and 2n using a Sieve of Eratosthenes. For each prime q between n and 2n, 
determine S{q), if it is less than g^'^^, by simply checking if /eg + 1 is prime for 
k = 1,2, ... , [g°-^^J. If we find a prime p = S{q) < g^-^^, add p to V. This 
is repeated until V contains (3i + (32 + i primes. If we are unable to find this 
number of primes, we have underestimated fi (since Theorem 4.2 guarantees 
their existence), so we double and restart the process. Obviously in practice 
we would not redo primality tests already performed for smaller n, so really no 
work need be wasted. 

Theorem 4.3. For log2 Ci < (3i, log2 C2 < (32, i, and n as in Theorem 4.2, we 
can generate (3i + (32+i elements ofV with 0{{(3i + (32+ey- \og^^°^^'^ {(3i + (32+i)) 
bit operations. At least C. of the primes in V will be useful. 

Proof. The method and correctness follows from the above discussion. The 
Sieve of Eratosthenes can be run with O (n log log log n) bit operations (see 
Knuth (1981), Section 4.5.4), and returns 0(n/logn) primes g between n and 
2n. Each primality test of /eg + 1 can be done with (logn)^"'""'^^) bit operations 
(Lenstra & Pomcrancc 2005), so the total cost is 0(n2(log n)^+°(^)) bit oper- 
ations. Since n G 0((/3i + /32 + ^) ■ log(/3i + (32 + i)) the stated complexity 
follows. □ 

The analysis of our methods will be significantly improved when more is 
discovered about the behavior of the least prime congruent to one modulo 
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a given prime, which we have denoted S{q). An asymptotic lower bound of 
S{q) E fi(glog q) is conjectured in Granville & Pomerance (1990), and we 
have employed the upper bound from Mikawa (2001) of S{q) G 0{q^'^^) (with 
exceptions). From our own brief computational search we have evidence that 
the conjectured lower bound may well be an upper bound: for all primes q < 
2^^, S{q) < 2gln^ q. If something similar could be proven to hold asymptotically 
(even with some exceptions), the complexity results of this and the next section 
would be improved significantly. In any case, the actual cost of the algorithms 
discussed will be a reflection of the true behavior of S{q), even before it is 
completely understood by us. 

Even more improvements might be possible if this rather complicated con- 
struction is abandoned altogether, as useful primes would naturally seem to be 
relatively plentiful. In particular, one would expect that if we randomly choose 
primes p directly from a set which has, say, 4(/3i + /52 + ^) primes, we might 
expect that the probability that p\Ci or {p — 1)\ C2 to less than, say, 1/4. 
Proving this directly appears to be difficult. Perhaps most germane results 
to this are lower bounds on the Carmichael Lambda function (which for the 
product of distinct primes pi, . . . , pm is the LCM of pi — 1, . . . , pm — I) , which 
are too weak for our purposes. See Erdos et al. (1991). 

5. Complexity analysis 

We are now ready to give a formal complexity analysis for the algorithms 
presented in Section 2 and Section 3. For all algorithms, the complexity is 
polynomial in the four bounds Ba, Bt, Bh, and B^ defined in Section 2.2, 
and since these are each bounded above by size(/), our algorithms will have 
polynomial complexity in the size of the output if these bounds are sufficiently 
tight. 

5.1. Complexity of Sparsest Shift Computation. Algorithm 2.3 gives 
our algorithm to compute the sparsest shift a of an unknown polynomial / G 
Q[x] given bounds Ba, Bt, Bh, and Bj^ and an oracle for choosing primes. 
The details of this oracle are given in Section 4. 

To choose primes, we set i = 2Ba + 1, and f3i = 2Bh and (32 = Bn{3Bt — 1) 
(according to Corollary 2.8). For the sake of notational brevity, define -Bs = 
Ba + Bh + BnBt so that Pi + + i e 0{Bj^). 

Theorem 5.1. Suppose f E Q[x] is an unknown polynomial given by a black 
box, with bounds Ba,Bt,Bh, and Bj^ given as above. If degf > 2Bt, then 
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the sparsest shift a G Q of f can be computed deterministically using 

bit operations, plus 0{KfB^^^ log^'*^ -BsM(log By;)) bit operations for the black- 
box evaluations. 

Proof. Algorithm 2.3 will always terminate (by satisfying the conditions of 
Step 2) after 2Ba + 1 good primes have been produced by the oracle. 

Using the oracle to choose primes, and because /3i + (32 + ^ ^ 
0{B§. log^"*""*-^-* By) bit operations are used to compute all the primes on Step 3, 
by Theorem 4.3. And by Theorem 4.2, each chosen p is bounded by 

All black-box evaluations are performed on Step 4; there are p evaluations 
at each iteration, and 0{By) iterations, for a total cost of O^KfBYp- M(logp)) 
bit operations. The stated complexity bound follows from the size of each 
prime p. 

Steps 10-14 are never executed when deg / > 2Bt. Step 15 is only executed 
once and never dominates the complexity. 

Dense polynomial interpolation over Zp is performed at most 0{By) times 
on Step 5 and 0{p) times at each of 0{Ba) iterations through Step 7. Since 
p ^ By, the latter step dominates. Using asymptotically fast methods, each 
interpolation of f^'P\x + 7) uses 0(M(p)logp) field operations in Zp, each of 
which costs 0(M(logp)) bit operations. This gives a total cost over all iterations 
of 0(-Byip Tog M(p) ■ M(logp)) (a slight abuse of notation here since the value 
of p varies). Again, using the fact that p G 0{B^^^ log^ *^^ -Be) gives the stated 
result. □ 

To simplify the discussion somewhat, consider the case that we have only a 
single bound on the size of the output polynomial, say Bf > size(/). By setting 
each of B^, Bh, and Bj^i equal to Bf, and by using the multiplication algorithm 
from Cantor & Kaltofen (1991), we obtain the following comprehensive result: 

Corollary 5.2. The sparsest shift a of an unknown polynomial f G Q[x], 
whose shifted-lacunary size is bounded by Bf, can be computed using 

O {Bf^ ■ log^-^8 Bf ■ (loglog Bf)'^ ■ logloglog Bf) 
bit operations, plus 

O {KfBy^ ■ log2-®9 Bf ■ loglog Bf ■ logloglog Bf) 
bit operations for the black-box evaluations. 
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Proof. The stated complexities follow directly from Theorem 5.1 above, 
using the fact that M(n) G O log loglog n) and G 0{B'j). Using the 
single bound Bf, we see that these costs always dominate the cost of Steps 10- 
14 given in Theorem 2.5, and so we have the stated general result. □ 

In fact, if we have no bounds at all a priori, we could start by setting Bf 
to some small value (perhaps dependent on the size of the black box or Kf), 
running Algorithm 2.3, then doubling Bf and running the algorithm again, 
and so forth until the same polynomial / is computed in successive iterations. 
This can then be tested on random evaluations. Such an approach yields an 
output-sensitive polynomial-time algorithm which should be correct on most 
input, though it could certainly be fooled into early termination. 

This is a significant improvement over the algorithms from our original pa- 
per (Giesbrecht & Roche 2007), which had a dominating factor of Bj^ in the 
deterministic complexity. Also — and somewhat surprisingly — our algorithm 
is competitive even with the best-known sparsest shift algorithms which require 
a (dense) f E Q[x] to be given explicitly as input. By carefully constructing the 
modular black box from a given / e Q.[x], and being sure to set Bt < (deg /)/2, 
we can derive from Algorithm 2.3 a deterministic sparsest-shift algorithm with 
bit complexity close to the fastest algorithms in Giesbrecht et al. (2003); the de- 
pendence on degree n and sparsity t will be somewhat less, but the dependence 
on the size of the coefficients log ||/|| is greater. 

To understand the limits of our computational techniques (as opposed to 
our current understanding of the least prime in arithmetic progressions) we 
consider the cost of our algorithms under the optimistic assumption that S{q) G 
0{q\r? q), possibly with a small number of exceptions. In this case the sparsest 
shift a of an unknown polynomial / G Q[x], whose shifted-lacunary size is 
bounded hj Bf, can be computed using 

O {B) ■ log^5; ■ {\og\og Bff ■ logloglogi?^) 

bit operations. As noted in the previous section, we have verified computation- 
ally that S{q) < 2q\v? g for g < 2'^^. This would suggest the above complexity 
for all sparsest-shift interpolation problems that we would expect to encounter. 

5.2. Complexity of Interpolation. The complexity analysis of the sparse 
interpolation algorithm given in Algorithm 3.5 will be quite similar to that 
of the sparsest shift algorithm above. Here, we need ^ = max{2BH + 1, -Bat} 
good primes to satisfy the conditions of Step 2, and from Theorem 3.7, we set 
/3i = 2BhBt and (32 = \B]^Bt{Bt — 1). Hence for this subsection we set 
Bs = Bt{Bh + BnBt) so that [3i+P2 + le 0{Bs). 
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Theorem 5.3. Suppose / G Q[x] is an unknown polynomial given by a mod- 
ular black box, with bounds Bt, Bh, Bn, and Bs given as above. The sparse 
representation of f as in (1.3) can be computed with 

O (^B 3 log Bs ■ MiB^-^' log'"'' Bs) ■ M (logics) 
+ BI,-M{{Bn + log Bt) \og{BN + log Bt))) 

bit operations, plus 0{KfB's^^ log^'*^^ i?5M(log Bs)) bit operations for the black- 
box evaluations. 

Proof. As in the sparsest-shift computation, the cost of choosing primes 
in Step 3 requires 0(i?| log^"*""*-^-* Bs) bit operations, and each chosen prime pk 
satisfies pk e 0{B\:^'nog^-^^ Bs). The total cost over all iterations of Step 4 
is also similar to before, 0{Bs ■ M{pk)logpk ■ M(logpA,.)) bit operations, plus 
0{Kf BspM [log pk)) for the black-box evaluations. 

We can compute each gf*^^'^ in Step 10 using 0(M(t) logt) ring operations 
modulo Pi — 1. Note that k G which is 0{Bh + B^), so the total cost 

in bit operations for all iterations of this step is 0{{Bh + Bjq) logi?T ■ M(i?r) ■ 
UilogBs)). 

Step 11 performs t Chinese Remainderings each of k modular images, and 
the size of each resulting integer is bounded by 2^'^ , for a cost of 0{Bt log B^ ■ 
M{Bj^)) bit operations. 

To factor g in Step 12, we can use Algorithm 14.17 of von zur Gathen & 
Gerhard (2003), which has a total cost in bit operations of 

O {Bl ■U{Bn + log Bt) + Bn{Bn + log Bt) ■ M((5jv + log Bt) log(5^ + log^r))) 

because the degree of g is t, g has t distinct roots, and each coefficient is 
bounded by 2^^ . 

In Step 15, we must first compute the modular image of mod pj — 1 and 
then look through all t exponents of f^^^^ to find a match. This is repeated tk 
times. We can use fast modular reduction to compute all the images of each 
Ci using 0{M{Biy) log -Bat) bit operations, so the total cost is 0{Bt{BhBt + 
Bj^Bt + M{Bn) log -Bat)) bit operations. 

Finally, we perform Chinese remaindering and rational reconstruction of 
t + 1 rational numbers, each of whose size is bounded by Bh, for a total cost 
oiO{BT-U{BH)logBH). 

Therefore we see that the complexity is dominated either by the dense 
interpolation in Step 4 or the root-finding algorithm in Step 12, depending 
essentially on whether B^ dominates the other bounds. □ 
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Once again, by having only a single bound on the size of the output, the 
complexity measures are greatly simplified. 

Corollary 5.4. Given a modular black box for an unknown polynomial 
/ e Q[x] and a bound Bf on the size of its lacunary represenation, that repre- 
sentation can be interpolated using 

O {BY'' log^-«9 Bj{\og\ogBf f logloglog 5^) 

bit operations, plus 

O {KfBf^'^ log^-®^ Bf loglog Bf logloglog Bf) 

bit operations for the black-box evaluations. 

Similar improvements to those discussed at the end of Section Section 5.1 
can be obtained under stronger (but unproven) number theoretic assumptions. 

6. Conclusions and Future Work 

Here we provide the first algorithm to interpolate an unknown univariate ratio- 
nal polynomial into the sparsest shifted power basis in time polynomial in the 
size of the output. The main tool we have introduced is mapping down modulo 
small primes where the sparse shift is also mapped nicely. This technique could 
be useful for other problems involving lacunary polynomials as well, although 
it is not clear how it would apply in finite domains where there is no notion of 
"size". 

There are many further avenues to consider, the first of which might be 
multivariate polynomials with a shift in each variable (see, e.g., Grigoricv & 
Lakshman (2000)). It would be easy to adapt our algorithms to this case 
provided that the degree in each variable is more than twice the sparsity (this 
is called a "very sparse" shift in Gicsbrccht et al. (2003)). Finding multivariate 
shifts in the general case seems more difficult. Even more challenging would be 
allowing multiple shifts, for one or more variables — for example, finding sparse 
gi, . . . ,gk G Q[x] and shifts ai, . . . , G Q such that the unknown polynomial 
/(x) equals gi{x — ai) + • • • + gk{x — a^)- The most general problem of this 
type, which we are very far from solving, might be to compute a minimal-length 
formula or minimal-size algebraic circuit for an unknown function. We hope 
that the small step taken here might provide some insight towards this ultimate 
goal. 
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